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ABSTRACT 

We present a general recipe for constructing A'-body realizations of galaxies comprised of 
near-spherical and disc components. First, an exact spherical distribution function for the 
spheroids (halo & bulge) is determined, such that it is in equilibrium with the gravitational 
monopole of the disc components. Second, an A'-body realisation of this model is adapted to 
the full disc potential by growing the latter adiabatically from its monopole. Finally, the disc 
is sampled with particles drawn from an appropriate distribution function, avoiding local- 
Maxwellian approximations. We performed test simulations and find that the halo and bulge 
radial density profile very closely match their target model, while they become slightly oblate 
due to the added disc gravity. Our findings suggest that vertical thickening of the initially 
thin disc is caused predominantly by spiral and bar instabilities, which also result in a radial 
re-distribution of matter, rather than scattering off interloping massive halo particles. 
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1 INTRODUCTION 

Generating an equilibrium A'-body representation of a multi- 
component galaxy is of importance for a number of applications, 
for example the study of bars (e.g. Debattista & Sellwood 2000; 
Athanassoula 2002), warps (e.g. Ideta et al. 2000) and galaxy merg- 
ers, both minor (e.g. Mihos et al. 1995; Walker et al. 1996) and 
major (e.g. Heyl, Hernquist, & Spergel 1996; Naab, Burkert, & 
Hernquist 1999). Yet, constructing realistic equilibrium models is 
considerably difficult and in fact impossible rigorously, i.e. requires 
some sort of approximation. 

Barnes (1988) introduced a "rather ad-hoc" method for con- 
structing A'-body models for a multi-component galaxy. He begins 
by constructing separate spherical equilibrium A'-body King (1966) 
models for the halo and bulge, which are then superposed and al- 
lowed to relax into a new equilibrium over several dynamical times. 
Next, the gravitational potential of the disc is slowly imposed, caus- 
ing both a flattening and a radial contraction of halo and bulge. Fi- 
nally the disc component is populated with particles drawn from 
a Maxwellian distribution with the local mean and dispersion in 
agreement with the Jeans equations. 

Hernquist (1993) introduced another method, which makes it 
possible to specify the density and velocity profiles of the various 
components in a straightforward way. In this method, the positions 
of all the particles in all components are determined first. This can 
be carried out easily since the desired density profiles are known. 
The Jeans equations are then used to find the velocity dispersions of 
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the various components, with those for the halo being found under 
the assumption that the potentials of disc and bulge can be approx- 
imated by their spherical average. Finally, individual velocities for 
all components are drawn from a locally Maxwellian distribution 
with the appropriate mean and dispersion velocities. 

Boily, Kroupa, & Penarrubia-Garrido (2001) extended this ap- 
proach to include a non-spherical halo, but maintained the approxi- 
mation to a Maxwellian distribution for the particle velocities. This 
method is not rigorous and Hernquist (1993) himself suggested that 
"In the future, it will likely be necessary to refine the basic approach 
as computer hardware and software permit simulations with parti- 
cle numbers significantly in excess of those discussed here." (he 
used N = 49 152 in his empirical tests). 

Kazantzidis, Magorrian, & Moore (2004b) compared the 
properties of A'-body haloes based on the Maxwellian approxima- 
tion with those derived from an exact distribution function found 
through Eddington (1916) inversion. As their study clearly demon- 
strates, the Maxwellian approximation results in non-equilibrium 
effects and is inappropriate, for example, when modelling the tidal 
stripping of substructure in CDM haloes. Unfortunately, Edding- 
ton inversion only works for spherically symmetric equilibria and 
cannot be straightforwardly applied to the problem of finding the 
self-consistent distribution function of a multi-component system. 

Therefore, Kuijken & Dubinski (1995) followed an alternative 
method, later expanded by Widrow & Dubinski (2005). They start 
from an analytical ansatz for the equilibrium distribution function 
of each component in isolation, written in terms of energy and an- 
gular momentum. In the case of the disc, the distribution function is 
also a function of the "vertical energy" E z = \y\ + <J?(K, z)-<J?(ft, 0), 
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which for orbits close to the disc plane is approximately conserved. 
These distribution functions are found for predetermined isolated 
density distributions. In the case of Kuijken & Dubinski (1995) 
these were a King (1966) model bulge, a lowered Evans model 
(Kuijken & Dubinski 1994) halo and an exponential disc. The com- 
bined model is then found by using these distribution functions in 
the combined gravitational potential of the entire model. The major 
problem with this approach is summed up by Widrow & Dubin- 
ski (2005) when they state that the combined model constructed in 
this way "may bare [sic] little resemblance to the corresponding 
isolated components, a situation which is cumbersome for model 
building". This approach is also rather inflexible. 

In this study, we therefore return to using Eddington inversion 
(in fact a generalisation thereof) for the near-spherical components 
but an analytic distribution function for the disc. Our method is 
detailed is Section 2, while Section 3 presents extended tests of 
the robustness and stability of the constructed A'-body equilibria. 
Finally, Section 4 sums up and concludes. 



2 THE NEW METHOD 

Our method bears some relation to that of Barnes (1988), in that the 
non-spherically symmetric component (the disc), is grown adiabat- 
ically within the A'-body halo. The bulge is assumed to be spher- 
ically symmetric, though it would be a relatively straightforward 
step to include a non-spherical bulge in a similar way to the disc. 
The process has three stages: 

(i) Creating an equilibrium initial A'-body representation of the 
spherically symmetric components of the galaxy (halo and bulge) 
in the presence of an external potential field corresponding to the 
monopole (spherical average) of the desired disc potential. 

(ii) Evolving the A'-body system by growing the non-monopole 
components of the disc potential adiabatically in a A'-body simula- 
tion, in order to allow halo and bulge particles time to relax into the 
new potential. 

(iii) Replacing the external potential field representing the disc 
component with an A'-body representation. 



2.1 Creating the initial A'-body halo and bulge 

We create an A'-body realization of the spherically symmetric com- 
ponents of the system using Cuddeford's (1991) method, an exten- 
sion of the model by (Osipkov 1979) and (Merritt 1985), which 
in turn extended Eddington's (1916) original inversion technique 
to obtain the distribution function for a spherical system with an 
isotropic velocity distribution. 

Cuddeford considered distribution functions of the form 



f(6,L) = L 2 "f (Q). 



(1) 



Here, Q = £ - l} j2r 2 a as for Osipkov-Merritt models (£ = *F - 
|v 2 , ¥ denoting the negative of the gravitational potential) with 
anisotropy radius r a . The parameter a is constrained to be greater 
than -1. As Cuddeford showed, f{Q) is related to the density by 
an Abel integral equation, which can (under the assumption that 
f(Q < 0) = 0) be inverted to yield 
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where n is defined as the largest integer equal to or less than a + ■ 
and where the "reduced density" is given as 
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This distribution function produces a spherically symmetric 
system with a velocity distribution such that 
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In the case where a = this model reduces to an Osipkov-Merritt 
model. In the case where r a — » 00, the anisotropy of the halo is the 
same at all radii, /? = -a. 

This approach has the advantage that the distribution function 
is exact for a spherically symmetric system, and thus remains in 
equilibrium, maintaining its original density profile. At no point is 
the assumption of a Maxwellian velocity distribution made. The 
only restrictions to this method are that *F = Y(r) is a monotonic 
function of radius (only); that p = p(r) and that the solution to 
equation (2) must be physical, i.e. fo(Q) > 0. This allows us to use 
a wide range of different halo (or bulge) density profiles. 

It is extremely useful that the derivation of equation (2) does 
not make the assumption that the potential of the system is that due 
to the density profile through the Poisson equation. This means that 
it is straightforward to generalise this approach to find the distri- 
bution function of a spherically symmetric component of a larger, 
spherically symmetric, system. The term p re <i in equation (3) is re- 
placed by the reduced density p red , the reduced density of the com- 
ponent ('; the term *P in equation (3) always refers to the potential 
of the entire system. 

In the model, the presence of the disc component is taken into 
account when calculating the distribution function of the halo and 
bulge. It is impossible to do this exactly using this method, as the 
disc is not spherically symmetric. The best approximation available 
in this case is to take the spherical average of the disc potential. 
Then the distribution functions of the bulge and halo are found in 
the total potential of the halo, bulge and fictitious spherically aver- 
aged disc. This prevents the radial contraction seen in the models 
constructed by Barnes (1988). 

Finally, we shall briefly describe how we choose initial con- 
ditions from a model with distribution function of the form (1). 
First, we draw an initial position from the density in the usual 
way by inverting the cumulative mass profile for the radius. Sec- 
ondly, we sample velocities by introducing pseudo-elliptical co- 
ordinates (u, 77) in velocity space such that v r = ucosr] and v, = 
w(l + r 2 /rl)~ [/2 sin;/. Then we pick values for rj and u at random 
from the distributions p(r]) = sin 1+2a 77 and p(u) = u 2+2a f Q C¥- ju 2 ), 
using the rejection method for the latter. 



2.2 Evolving the halo and bulge to adapt to the disc 

The second stage of creating the model is to evolve A'-body model 
of halo and bulge so that it is in equilibrium with the full disc po- 
tential, rather than with its spherical average. To this end we use 
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the iV-body code gyrfalcON, which is based on Dehnen's (2000; 
2002) force solver falcON — though the axisymmetry of the sys- 
tem means that a code based upon an expansion in spherical har- 
monics about the origin (such as the so-called self-consistent field 
codes, e.g. Hernquist & Ostriker 1992) would be well suited to this 
purpose. 

The full potential of the disc, Odi sc (*) is grown from its spher- 
ical average O disc (r) according to the formula 



<b(x, t) = O disc , (r) + A{t) [O dlsc (x) - OdiscoM] , 



(7) 



where A{t) is a growth factor that goes from at t = to 1 smoothly 
in a time-scale, f grow , far longer than the dynamical time in the re- 
gion of the disc. The halo and bulge are then allowed further time 
to relax completely under the influence of the disc potential. 

This process causes changes to the density distribution of the 
halo. The spherically averaged density distribution is largely un- 
changed, but the halo and bulge are somewhat flattened, and the 
iso-density contours of both the halo and bulge become oblate. The 
degree of flattening is dependent on the details of the various com- 
ponents, a typical example is given in Section 3. 

This stage is by some distance the slowest in the creating of 
the initial conditions. However it is still takes a relatively small 
fraction of the total CPU time of almost any scientifically interest- 
ing simulation. It should also be noted that up until this point the 
velocity distribution of the disc has not been a factor in the calcula- 
tions, so the halo and bulge created by this process can be re-used 
in simulations with identical disc density profiles, but different disc 
kinematics. 



2.3 Populating the disc 

We start from the assumption that the motion of disc particles 
decouples into its components in the plane of the disc, and per- 
pendicular to it. That is, we assume that the planar and vertical 
components of the total energy, E\\ = ^(v 2 R + v?) + OCR, 0) and 
E ± = |v? + ®(K, z) - $>(R, 0), are both separately conserved. This 
assumption is usually excellent for orbits near the disc, for which 
always \z\ <k R. The potential O here is, of course, the total potential 
of the system, i.e. that of the disc model plus that of the evolved N- 
body system. The first is computed as in Dehnen & Binney (1998), 
while the latter is approximated using a potential expansion (so- 
called SCF) method whereby ensuring axial symmetry by taking 
only m = and even / terms. 

We make an ansatz for the disc distribution function, following 
the approach of Dehnen (1999b). As in that study, we first consider 
the simplest distribution function, i.e. that of a dynamically com- 
pletely cold disc, in which all particles are on circular orbits. For a 
disc with density p = I.(R) 6(z) (X being the surface density), this 
can be written as (e.g. Dehnen 1999b) 

Sl(R L )X(R L ) 

f coW (E h L z , z, v z ) = j— — <5(£|, - E c (R Lz )) S(z) 6(v z ), (8) 

where C1(R) and k(R) denote the angular and epicycle frequency 
at radius R and E C (R) the energy of the circular orbit through R, 
while R L _ is the radius of the circular orbit with angular momen- 
tum L z . In order to obtain a warm disc distribution function, one 
may replace the ^-functions in (8) by exponentials, resulting in a 
locally Maxwellian velocity distribution (e.g. Shu 1969; Kuijken & 
Dubinski 1995). 

However, as demonstrated by Dehnen (1999b), it is better to 



use the following alternative form for the cold-disc model before 
"warming" it up. 

/ cold (£|j,L z ,z,v z ) = 



m<(R E ..) 



■S(n(R Et )[L-L c (R Et )])6(z)6(v z ),(9) 



where L C (R) is the angular momentum of the circular orbit with 
radius R and R E the radius of the circular orbit with energy E\\. 
Again, the warm disc distribution function is obtained by replacing 
the (5-functions with exponentials (noting that z = 0, v, = means 
that E ± =0), giving 

Q(R E )t(R E ) 1 
U x (E h E ± ,L z ) = ——exp 
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(10) 



Here, E(/f) and a R (R) are sought such that the true surface-density 
and radial velocity-dispersion profiles of the A'-body representation 
are those desired, to within an appropriate degree of accuracy, while 
the velocity dispersion in the vertical direction cr 2 z (R) = nGzd ^(R) 
with disc scale height parameter za- 

Dehnen (1999b) argued that equation (10) gives a more use- 
ful distribution function than that of Shu (1969); he pointed that 
the warming of a disc can be described as an exponential in the ra- 
dial action J R . The radial action is more closely approximated by 
Cl(R Et )[L c (R Et )-\L z \y K (R Et ) than by [E\\- E C (R L: )M k(R L: ) (Dehnen 
1999a). More practically, this distribution function has the advan- 
tage that the value of R E is generally a far better approximation to 
the mean radius of an orbit than R L , , which ensures that ~L(R) and 
cr R (R) closely resemble £(/?) and cr R (R), respectively. This choice 
of distribution function also extends to negative L z , unlike that of 
Shu, allowing for a tail of counter-rotating stars. 

When sampling from this distribution function, we attempt 
to minimise noise in the particle distribution by sampling points 
more regularly than random in phase space. We accomplish this 
by sampling orbits from the density distribution Y.(R), then placing 
1 <k N sam <<c JV disc particles at points on each orbit. Finally, we iter- 
atively adapt i.(R) and d- R (R) so that the actual surface density and 
velocity dispersion profile match the target profiles as closely as 
possible. The procedure for sampling the planar part of the phase- 
space positions is very similar to that proposed in Dehnen (1999b) 
and its details are as follows. 

(i) Draw a radius from a thin disc with surface density f,(R) (ini- 
tially S = S) and set E l{ to the energy of the circular orbit at R. 

(ii) Draw a number £ 6 (0, 1) and determine 

L z = L c (R) + ln^di(R)/Q(R). (11) 
If L z t [-L C (R\ L C (R)] go back to step (i). 

(iii) Integrate the orbit with these values of £y and L z for one 
radial period T R , find the radial frequency oj r , and tabulate R(t) and 
R(t) for that orbit. Evaluate the correction factor g corr = k(R)/w r 
(see Dehnen 1999b). 

(iv) Given N olb ss N 1 ^, choose N sam to be either of the two in- 
tegers next to gcorr-Norb drawn with probabilities such that the mean 
equals g cm rN olb . Next sample A^ sam orbital phases f,- 6 (0, T R ) and az- 
imuth angles 0, 6 (0, 2n), and determine the corresponding phase- 
space points {R h <pj,Ri, <f>i) from the orbit. 

(v) Repeat steps (i) to (iv) until a total of A^ disc disc particles have 
been sampled. 
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Figure 1. Rotation curve for the test galaxy model. The solid line is the 
net rotation curve, also shown are the contributions from the disc (dotted), 
bulge (short-dashed), and halo (long-dashed). 



(vi) Finally, evaluate from the N Asc phase-space points the actual 
surface density and radial velocity dispersion, Z ou , and cr oulR , of the 
N-body model and adapt 



W) 

SoutW 



cr R (R) -> cr R (R) 



o-r(R) 

CToutJ>(R) 



(12) 



Then repeat the whole sampling procedure of steps (i) to (v) and 
iteratively adapt £(/?) and d- R (R) until no further improvement oc- 
curs. 

If we use quasi-random, rather than pseudo-random, numbers 
in this procedure, it closely resembles the "quiet start" method 
(e.g. Sellwood 1987). Once this is done we determine the verti- 
cal component of each disc particle's position and velocity. We 
assume that the local structure in the z-direction corresponds to 
that of an isothermal sheet with a constant vertical scale height, 
H (Spitzer 1942). This leads to a density throughout the disc 
p(R,z) oc sech 2 (z/z d ). v z is then drawn randomly from a normal 
distribution with a 1 , = 7iGI.(R E ) z<i- 



bulge mass M b = 0.2 and scale length r b = 0.2. Scaling these values 
to the Milky Way, taking R d = 3.5 kpc, M d +M b = 5x 10 10 M Q (as in 
Klypin et al. 2002) gives a time unit ^ 14 Myr, and thus a velocity 
unit of ~ 250 km s" 1 . 

For the halo we take a scale radius r h = 6, truncation radius 
r t = 60, and mass M b = 24, 79% of which is within the truncation 
radius. The rotation curve for this model (with velocities given in 
code units) is shown in Figure 1. 

In all tests in which they were populated, the halo had 
1200000 particles, the bulge 40000 and the disc 200000. This 
corresponds to each halo particle being 4 times more massive 
than a stellar particle. We use gravitational softening lengths of 
6 = 0.02 for disc and bulge particles, and 0.04 for halo particles. 
This choice ensures that the maximum force exerted by a single 
particle (oc nij/ef ) is the same for all particles. 

Simulations were performed with gyrfalcON with a minimum 
time-step of 2~ 7 ', and a block-step scheme with largest time step 2" 4 . 
Individual particle time-steps were adjusted such that on average 
the time-step 



t, = mm 



0.01 
IfliT' 



0.05 



(16) 



with O, and a, the gravitational potential and acceleration of the ith 
body in simulation units. 

We first tested that the spherically symmetric components re- 
main in equilibrium in the spherically-symmetrised potential. This 
was tested with various different choices for r a and a (equation (1)) 
in the halo, in order to ensure that Cuddeford inversion had been 
implemented correctly. The halo and bulge profiles remained con- 
sistent to within a few softening lengths of the centre, where soft- 
ening and two-body relaxation have some small effect. This testing 
was also used to assess whether the time integration parameters 
were appropriate. Since energy was typically conserved to within 
0.05% over 200 time units in these test simulations, we assume that 
the time integration is sufficiently accurate. For all further tests we 
only consider the isotropic case (a = 0, r a — > oo, so/3 = through- 
out the halo). 



3 TESTING THE NEW METHOD 

We test our approach using a model based loosely upon the Milky 
Way as modelled by Klypin, Zhao, & Somerville (2002). The disc 
is defined as having an exponential surface density profile with a 
vertical structure modelled by isothermal sheets, i.e. 



Pa(Rz) 



4^z d 



exp 



- — ) sech* £ 



(13) 



where M A is the total disc mass, R A is the disc scale radius, and Zd 
is a scale height (though it should be noted that this is not the e- 
folding height). The halo model is a truncated (Navarro et al. 1997, 
hereafter NFW) halo with initially spherical density distribution 

sech(r/r t ) 



Ph(r) = p c 



(r/r h )(l + r/r h y' 



(14) 



while the bulge is modelled with a Hernquist (1990) density profile 
M b r b 



Pb(r) = 



(15) 



2nr(r h + r) 3 

We choose units such that Newton's constant of gravity G = 1, 
Rd = I, and M d = 1. We take the disc scale height Zd = 0.1, the 



3.1 Testing the disc growth 

The next step was to ensure that the growth of the full disc poten- 
tial from its monopole occurs without significant effect upon the 
radial density profile of the halo or bulge, or upon their kinematics, 
and to quantify the effect upon the shape of the resultant density 
distribution. The disc potential is grown as per equation (7), with a 
total growth time t grow = 40. This is substantially longer than the 
dynamical time f dyn in the vicinity of the disc (for instance at R = 3, 



'dyn 



~6). 



Figure 2 shows the evolution of the density profiles of halo 
and bulge. Both are maintained to within a high degree of accuracy 
through the growth of the full potential, and then remain close to 
unchanged for a long period of time. The halo density is slightly 
raised in the inner ~ 0.02r h = 0.12 and the bulge density is slightly 
lowered in the inner ~ 0.3/1, = 0.06, in each case corresponding 
to only a few softening lengths. This is approximately the same as 
observed in simulations with no disc growth and is likely caused by 
mass segregation due to (artificial) two-body relaxation. 

To quantify the flattening of halo and bulge due to the growth 
of the full disc potential, we determine the axis ratios of the halo 
and bulge as a function of radius. To this end, we first estimate the 
local density at each particle's position by the nearest-neighbour 
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Figure 2. Spherically averaged density profile of the halo (left) and bulge (right), shown just before growth of the full disc potential (t = 0), and at 100 and 200 
time units. The disc component was not populated in this simulation. The density profiles remain more or less unchanged by the growth of the full potential, 
except at radii only slightly larger than the respective softening length (eh = O.OO6601 and e\, = O.ln,, indicated by the dotted vertical line). 




r A h r / r t, 

Figure 3. Minor and intermediate axis ratios (cja < b/a respectively), for both the halo (left) and bulge (right) in the same simulations as in Fig. 2. Particles 
are binned in density, and the values plotted here are the axis ratios of these bins against their median radius. Deviations from sphericity (cja = h/a = 1) at 
t = are entirely due to discreteness noise. 



method of Casertano & Hut (1985), using the 15 nearest neigh- 
bours (of the same component). Next, we bin the particles in esti- 
mated density and evaluate the axis ratios in each bin as the square 
roots of the ratios of the eigenvalues of the moment-of-inertia ten- 
sor 1 . Figure 3 plots the axis ratios such obtained for halo and bulge 
versus the radius (the median radius of the corresponding density 
bin). Obviously, there is some small effect on the bulge and halo 

1 In the literature one often finds axis ratios computed for bins in spher- 
ical radius (e.g. Kazantzidis et al. 2004a). This not only requires accurate 
knowledge of the centre position, but more importantly results in a sub- 
stantial bias (because of the usage of spherical symmetry). For a triaxial 
body, for example, the axis ratios at small radii are drastically overestimated 
(Athanassoula 2007). Binning in potential energy gives somewhat better but 
still unsatisfactory results (since the gravitational potential is less flattened 
than the density). 



shape from the growth of the disc. In both cases the shortest (c) 
axis is very close to the z-axis as expected. The degree of flattening 
is relatively modest, with c/a - 0.8 in the most flattened shells. In 
the case of the halo this is within the inner ~ 0.3r h = 1.8i?j. In the 
case of the bulge the most flattened shells are at ~ 3rb = 0.6Ra. 
The bulge is less flattened at smaller radii because the disc is of 
finite thickness (scale height z d = 0. IRd), which means that the full 
potential is less flattened - when compared to the spherical average 
- in the inner parts of the bulge than in the outer. The same effect 
would be observable in the halo with sufficient resolution. 

For both halo and bulge Figure 4 shows the evolution of the 
radial profile of the velocity anisotropy parameter /? as determined 
from spherical shells. Both components are initially isotropic and, 
to within numerical accuracy, remain so during and well after the 
growth of the full disc potential. 
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r/r h 

Figure 4. Anisotropy parameter fi, determined from spherical shells, plotted 
against radius for the halo and bulge in the same simulations as in Fig. 2. 
Initial conditions were defined as having ji = and variation from that at 
t = is due to discreteness noise. 



3.2 Testing the disc model 

Before testing the full compound galaxy, we test the disc model in 
isolation, i.e. we follow the evolution of the populated disc com- 
ponent in the static external potential 0(7?,z), the cylindrically- 
symmetrised potential in which the distribution function of the disc 
was constructed. The purpose of this test is to ensure that the ap- 
proximate nature of our disc distribution function (the assumption 
that the planar and vertical energies are separately conserved) does 
not compromise our approach, i.e. that the disc A'-body model is 
close to equilibrium. 

The main difficulty here is that we want the disc model to be 
stable to axisymmetric perturbations, but not to non-axisymmetric 
instabilities, such as spiral waves and bar modes, both of which one 
would expect to see. These instabilities cause redistribution of the 
mass of the disc in both the R and z directions (e.g. Hohl 1971; 
Athanassoula & Misiriotis 2002). In the tests presented here, we 
prevent non-axisymmetric modes from growing by a technique pi- 
oneered by Athanassoula (private communication): the azimuth of 
every disc particle is randomised after every block-step, thus de- 
stroying any coherent non-axisymmetric perturbation. 

The distribution function defined by equation (10) allows for 
many possible cr R (R). In practise our current implementation re- 
stricts the possibilities to either cr R oc expi-R/R^), or cr R is such 
that the Toomre (1964) stability parameter 

Q S ™ (17) 
^ 3.36 GS v ' 

is constant throughout the disc (this should not be confused with 
the Osipkov-Merritt Q in equation (1)). A stellar disc is known to 
be unstable to axisymmetric waves if Toomre's Q < 1. 

We test two models, one with constant Q = 1.2, and one with 
cr R oc exp(-R/2Ra) (i.e. Ra- = Ra/2 so that cr R oc E), with the 
constant of proportionality defined such that Q(R = R^) = 1.2. The 
two give qualitatively similar results, so only the results from the 
latter distribution function are shown in Fig. 5 for simplicity. The 
surface density of the disc model is preserved to within the inner 
~ 0.1 -0.2R d , and out to beyond 6R d . The radial velocity dispersion 
cr R is preserved over a similar range. 

In the middle panel of Fig. 5, we plot the r.m.s. value of z as 
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Figure 5. Radial profiles of surface density (top), thickness (middle), 
and radial velocity dispersion (bottom) for an isolated disc model. The 
randomising-azimuth method has been used to prevent the growth of bar 
or spiral modes. 



a measure of the disc thickness. For a disc with p oc sech 2 (z/z d ) 
one expects z rms = kzaI Vl2 ss 0.907z,j, indicated by a dotted line. 
The disc thickness remains near constant in the range 0.3 < R < 4, 
with a slight warming which can reasonably be attributed to particle 
softening. The disc becomes somewhat flattened in its very inner 
(~ 03R d ) and outer (R > 4R d ) parts. This is presumably caused by 
too low initial z-velocities, which were assigned assuming that the 
vertical force is dominated by the local disc. In the inner parts of 
the disc, the bulge has a non-negligible contribution to the vertical 
force field, and in the outer parts the local disc is very tenuous, 
so the contribution of the halo (and the monopole of the disc) to 
the vertical force field is significant. It should also be noted that 
the approximation that the planar and vertical motion decouple is 
worst in the inner parts of the disk, which may well contribute to 
the flattening there. 



3.3 Testing the full model 

3.3.1 A constrained model 

Finally we test the complete compound model. First we wish to 
examine the behaviour of the system in the absence of bar or spi- 
ral instabilities. In order to do so, we utilise the same method of 
randomising the azimuth (of disc particles only) as in Section 3.2. 

However, using this method with a live halo and bulge is 
not straightforward, since a live A'-body simulation will often drift 
slightly from its original centre. When this happens the disc will 
be pulled in same the direction as the inner halo, so that randomis- 
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Figure 6. Spherically averaged density profile of the halo (left) and bulge (right) in a fully populated simulation with non-axisymmetric perturbations sup- 
pressed (see text). The profiles are shown at t = 0, when the full /V-body model was populated but before it has evolved; and at two times during the full 
rV-body simulation. Symmetry about the origin was enforced, to avoid numerical difficulties relating to off-centring. 
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Figure 8. Anisotropy parameter /3, as determined from spherical shells, vs. 
radius for the halo and bulge in the same simulations as Figs. 6 & 7. 



ing azimuths around the origin would alter the model structure. In 
order to prevent such a drift, we employed another symmetrisation 
method: the distribution of particles in halo and bulge was kept re- 
flection symmetric about the origin at all times. To this end, we 
treated particles in bulge and halo as pairs and enforced, initially 
and after each time step (during both the disc growth and the fully 
populated simulation), that for each pair positions and velocities 
both add up to zero. 

In this test we set the velocity dispersion <r R of the disc such 
that Toomre's Q = 1.2 at all radii. The full disc potential was grown 
over a period ? grow = 40, as in Section 3 . 1 , and then kept in place for 
a time fh id = 20 to ensure that the halo and bulge had fully relaxed 
in its presence. Only then was the disc populated with particles. 
Simulations were then run to observe the evolution over another 
200 time-units. Times quoted in Figs. 6-9 take t = to be the time 
when the fully populated simulation begins (not the start of the 
growth of the full disc potential, as in Section 3.1). 

Figures 6 and 7 show that the properties of the halo and 



bulge are maintained to within a high degree of accuracy through- 
out the simulation. The halo density is slightly raised in the inner 
~ 0.02r h = 0.12 and the bulge density is slightly lowered in the 
inner ~ 0.3r b = 0.06, much as they were in simulations with an 
unpopulated disc (Fig. 2). Both components remain isotropic (with 
some noise, Fig. 8), and while they are non-spherical, they do not 
become significantly more or less aspherical over the course of the 
simulation (Fig. 7). 

Figure 9 (top panel) shows that the surface density of the disc 
is nearly unchanged at all radii throughout the entire simulation. 
Also the disc thickness remains near constant over the full radial 
range. The small deviations are similar to those observed for the 
disc in isolation (Fig. 5). In particular, the disc does not thicken 
up as one might expect from heating due to collisions with (more 
massive) interloping halo particles. This is probably a side effect of 
the randomisation of the azimuthal angle of the disc particles. This 
randomisation has the effect that disc particles are rarely close to 
any halo particle crossing the disc for more than one time step. 

The radial velocity dispersion cr R in the disc (Fig. 9, bottom 
panel) is nearly unchanged in the range 0.5 < R < 3. In the outer 
parts of the disc there is some variation, but no consistent warming 
or cooling. In the inner parts of the disc there is a clear increase in 
cr R . This warming clearly does not affect the density distribution. It 
can reasonably be attributed to the facts that the decoupling of the 
vertical and planar motion is a poor approximation in the inner disc, 
and that, while the choice of a distribution function with constant Q 
can be useful in that it avoids the inner parts of the disc being very 
hot, it does lead to the velocity dispersion being unrealistically low 
at the very inner radii. 

3.3.2 Unconstrained tests 

Finally we run simulations in which the galaxy model is uncon- 
strained from the moment the disc component is populated. This 
means that the disc component is free to develop bar modes and 
other instabilities. 

In Fig. 10 we plot the surface density profiles of discs with var- 
ious initial cr R profiles. In each case cr R was defined such that Q was 
constant across the disc. Simulations are shown with Q = 1.2, 2, 3, 
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Figure 7. Minor and intermediate axis ratios (c/a < bja respectively), for both the halo {left) and bulge (right) in the same simulations as in Fig. 6. Particles 
are binned in density, and the values plotted here are the axis ratios of these bins against their median radius. 
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Figure 13. Evolution of the disk thickness in four different simulations: constrained (by randomising the azimuth) and unconstrained, with rigid or live halo. 
In all cases, the initial velocity dispersion of the disc was set such that Q = 1.2 everywhere. The lower-left panel corresponds to the middle panel of Fig. 9 (but 
note the different scales). 



and 4. For hotter discs (larger Q), it becomes increasingly difficult 
to iterate the particle positions such that the surface density is close 
to that desired, using the method described in Section 2.3. This is 
why the initial surface density of the Q = 4 disc (and to a lesser 
extent the Q = 3 disc) is not an exponential in radius. It should 
also be noted that the motivation of warming an initially cold dis- 
tribution function is increasingly invalid as the disc becomes hotter. 
Moreover, such hot axisymmetric discs are less realistic, given the 
observed fraction of barred galaxies (which cannot form from such 
hot discs). 

It has long been recognised that numerical simulations of cool 
discs are much more susceptible to bar formation than those of 
warm discs (e.g. Hohl 1971). This is what we observe in our mod- 
els. The scatter plot of particle positions for the Q = 1.2 case, Fig- 
ure 11, shows spiral structure forming within 50 time units, and the 
formation of a strong bar within 100. This has the effect of substan- 
tially altering the surface density profile, dramatically increasing 
the disc's surface density in the inner ~ 0.3 scale lengths (as seen 



in Fig. 10). The Q = 3 disc shows no signs of developing spiral 
structure within 200 time units (Fig. 12), and the density profile 
remains nearly unchanged throughout the simulation. 



3.3.3 On the origin of disc thickening 

It is instructive to compare the disc thickening in the four different 
cases we have looked at in this study (populated/rigid halo; ran- 
domised azimuth/unconstrained model), see Figure 13. Simulations 
which are constrained so that non-axisymmetric modes are sup- 
pressed show hardly any disc thickening, strongly suggesting that 
the dominant cause of disc thickening in our models is the action of 
non-axisymmetric modes within the disc, rather than the impact of 
halo particles crossing the disc. Athanassoula (2002) showed that a 
populated halo can stimulate bar growth by absorbing angular mo- 
mentum (which a rigid halo cannot). This is a likely explanation 
for the greater increase in disc thickness seen in the unconstrained 
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Figure 9. Radial profiles of surface density (top), thickness (middle), and 
radial velocity dispersion (bottom) of the disc component in a full W-body 
model, whose halo and bulge properties are shown in Figs. 6 to 8. The 
randomising-azimuth method has been used to prevent the growth of bar or 
spiral modes (see text). 




Figure 10. Surface density profiles for the discs in unconstrained simu- 
lations with four different values of Toomre's Q initially. Results for the 
different simulations are offset by a constant. 
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Figure 11. Scatter plot of disc particle positions in the x-y plane at various 
times in an unconstrained simulation with initial Q = 1.2. Only 1 in every 
20 particles is plotted, in the interests of clarity. 
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Figure 12. Scatter plot of disc particle positions in the x-y plane at various 
times in an unconstrained simulation with initial Q = 3. As in figure 11, 
only 1 in every 20 particles is plotted. 



simulation with live halo compared to that with a rigid halo (right 
panels in Fig. 13). 

Possibly the randomisation of the disk particles' azimuths is 
also a direct cause of a reduction in disc thickening, because of he 
rapid jumps it causes in the disc particles' positions (rather than 
just an indirect cause, through preventing non-axisymmetric insta- 
bilities). Why this would be the case, however, is unclear. 
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4 CONCLUSIONS 

We have described and tested a new method for constructing an 
equilibrium A'-body representation of a galaxy with halo, disc and 
bulge components. 

We have used a distribution function for the halo and bulge 
based upon Cuddeford (1991) inversion, while the distribution 
function of the disc is based on the work of Dehnen (1999b). One 
advantage over previous methods is that our method avoids the use 
of a Maxwellian approximation, and produces models which tend 
to stay very close to their original states, though non-axisymmetric 
instabilities develop in the disc if it is unconstrained and reasonably 
dynamically cold. 

Another advantage of our method is that the density distribu- 
tions of the various components are straightforward to prescribe. 
That is, the way the models are constructed guarantees that the 
equilibrium A'-body model has bulge, halo, and disc components 
the density profiles (as well as the velocity-anisotropy profiles) 
of which follow the target models very closely. Having said that, 
we have little control with our method over the degree of non- 
sphericity of halo and bulge, introduced by the disc's gravitational 
potential. We find that typically bulge and halo are mildly oblate 
(axis ratio ~ 0.8) in the region dominated by the disc. 

The choice of disc distribution function, while physically mo- 
tivated, and clearly in equilibrium (though potentially unstable to 
non-axisymmetric modes), causes problems when creating a warm 
disc (Toomre's Q > 4) since it becomes increasingly difficult to 
tailor the distribution function to the desired density and velocity 
dispersion profiles. However, we are not aware of any other method 
to create such warm A'-body disc with a truly exponential surface 
density profile. 

We performed several tests to validate that the A'-body model 
created by our method meets our expectations. These tests suggest 
that the growth of the disk thickness is predominantly driven by 
non-axisymmetric modes in the disc itself rather than direct inter- 
actions with halo particles. 

The computer programs generated in the course of this project 
will be publicly available under the NEMO computer package 
(www.astro.umd.edu/nemo). 
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